#Tomz data analysis
set.seed(5123)
working.dir <- "/Users/Allan/Dropbox/!!Papers/1Reputation/Honor-Culture-War/2012-05/Replication/Tomz/"
setwd(working.dir)

require(foreign)
require(ri)
data <- read.csv("130108_tomz.csv")
data <- data[complete.cases(data),]
#dropped one observation

#e1_cond==1 is no threat, e1_cond==2 is a threat. 

#subsetting on while males
r <- ppgender=="male"  & ppeth=="white, non-hispanic" & (e1_cond=="p1" | e1_cond=="p2")

data <- data[r,]

attach(data)
names(data)

r <- ppgender=="male"  & ppeth=="white, non-hispanic" & (e1_cond=="p1" | e1_cond=="p2")

length(South[e1_cond=="p1" & r])
length(South[e1_cond=="p2" & r])

s.nt <- dis[South==1 & e1_cond=="p1" & r]
ns.nt <- dis[South==0 & e1_cond=="p1" & r]
s.t <- dis[South==1 & e1_cond=="p2" & r]
ns.t <- dis[South==0 & e1_cond=="p2" & r]

testd0 <- (mean(s.nt)-mean(ns.nt))
#testd0 <- (mean(dis[South==1 & e1_cond=="p1" & r])-mean(dis[South==0 & e1_cond=="p1" & r]))
testd1 <- (mean(s.t)-mean(ns.t))
#testd1 <- (mean(dis[South==1 & e1_cond=="p2" & r])-mean(dis[South==0 & e1_cond=="p2" & r]))
testdd <- (mean(s.t)-mean(ns.t)) - (mean(s.nt)-mean(ns.nt))
#testdd <- (mean(dis[South==1 & e1_cond=="p2" & r])-mean(dis[South==0 & e1_cond=="p2" & r])) -
  (mean(dis[South==1 & e1_cond=="p1" & r])-mean(dis[South==0 & e1_cond=="p1" & r]))
logodds <- function(x) {
  result<-log(mean(x)/(1-mean(x)))
  return(result)
}
testdlor <- (logodds(s.t)-logodds(ns.t))- (logodds(s.nt)-logodds(ns.nt))
  
perms <- 10000
pS <- matrix(NA, nrow=length(South), ncol=perms)
i <- 1
while(i < perms+1){
  pS[,i] <- sample(South, length(South), replace=FALSE)
  i <- i+1
}
origpS <- pS
pS <- unique(origpS, MARGIN=2)

length(South[e1_cond=="p1" & r])
testsd0 <- rep(NA, length(perms))
testsd1 <- rep(NA, length(perms))
testsdd <- rep(NA, length(perms))
testsdlor <- rep(NA, length(perms))

i <- 1
while(i < perms+1){
  s.nt <- dis[pS[,i]==1 & e1_cond=="p1" & r]
  ns.nt <- dis[pS[,i]==0 & e1_cond=="p1" & r]
  s.t <- dis[pS[,i]==1 & e1_cond=="p2" & r]
  ns.t <- dis[pS[,i]==0 & e1_cond=="p2" & r]
  testsd0[i] <- (mean(s.nt)-mean(ns.nt))
  #testsd0[i] <- (mean(dis[pS[,i]==1 & e1_cond=="p1" & r])-mean(dis[pS[,i]==0 & e1_cond=="p1" & r]))
  testsd1[i] <- (mean(s.t)-mean(ns.t))
  #testsd1[i] <- (mean(dis[pS[,i]==1 & e1_cond=="p2" & r])-mean(dis[pS[,i]==0 & e1_cond=="p2" & r]))
  testsdd[i] <- (mean(s.t)-mean(ns.t)) - (mean(s.nt)-mean(ns.nt))
  #testsdd[i] <- (mean(dis[pS[,i]==1 & e1_cond=="p2" & r])-mean(dis[pS[,i]==0 & e1_cond=="p2" & r])) -
   # (mean(dis[pS[,i]==1 & e1_cond=="p1" & r])-mean(dis[pS[,i]==0 & e1_cond=="p1" & r]))
  testsdlor[i] <- (logodds(s.t)-logodds(ns.t))- (logodds(s.nt)-logodds(ns.nt))
  i <- i+1
}
testsd0 <- c(testd0, testsd0)
testsd1 <- c(testd1, testsd1)
testsdd <- c(testdd, testsdd)
testsdlor <- c(testdlor, testsdlor)

#Two sided test
if (rank(testsd0)[1]/(perms+1) <0.5){
pd0 <- (rank(testsd0)[1]/(perms+1))*2
}
if (rank(testsd0)[1]/(perms+1) >=0.5){
  pd0 <- (1-(rank(testsd0)[1]/(perms+1)))*2
}

pd0

pd1 <- 1 - rank(testsd1)[1]/(perms+1)
pd1

pdd <- 1 - rank(testsdd)[1]/(perms+1)
pdd

pdlor <- 1 - rank(testsdlor)[1]/(perms+1)
pdlor

